library(sf)
library(tidyverse)
#インタラクティブな地図作成
library(mapview)ある地物(地上にある物)や地点周辺の統計地図の作成,および統計量(合計値,平均値など)の計算.ここでは,多摩ニュータウン(ポイントデータ)周辺の人口(メッシュデータ)を可視化.また多摩ニュータウンから500mの円を描き,その円と重なるメッシュデータの人口数を計算.
知りたい統計量の例
行政区域
read_sf()):行政区域をTokyo_mapと名付ける.
filter()):Tama_shi_mapと名付ける.#東京の行政区域(令和4年)
Tokyo_map<-
read_sf("N03-20220101_13_GML/N03-22_13_220101.shp")
#多摩市周辺(多摩,町田,日野,稲城,八王子)
Tokyo_map %>%
filter(N03_007=="13224" | N03_007=="13209" |
N03_007=="13212" | N03_007=="13225" |
N03_007=="13201") ->
Tama_shi_map
#多摩市の行政区域の可視化
ggplot()+
geom_sf(data=Tama_shi_map)ニュータウン
read_sf()):行政区域をTokyo_newtownと名付ける.crs(座標参照系)を設定.
WGS84:1984年に定められた世界測地系.Tama_newtownと名付ける.#東京ニュータウン
Tokyo_newtown<-
read_sf("P26-13_13/P26-13_13/P26-13_13.shp", crs="WGS84")
#多摩ニュータウンの選抜
Tokyo_newtown %>%
filter(P26_005=="多摩ニュータウン") ->
Tama_newtown統計データ:500mメッシュ
市区町村別メッシュ・コード一覧(残念ながら3次メッシュの情報だけの模様)から多摩市のコードが5339からはじまることがわかる.
人口
read.table()を用いてデータを読込.Tokyo_popと名付ける.#人口
Tokyo_pop<-read.table("tblT001101H5339/tblT001101H5339.txt",
header=TRUE, sep=",")
#便宜上1行目削除(データ名が2行にわたるため)
Tokyo_pop %>%
slice(-1) ->
Tokyo_pop
#メッシュデータと結合するため,KEY_CODEのデータ型変換
Tokyo_pop %>%
mutate(KEY_CODE=as.character(KEY_CODE)) ->
Tokyo_pop境界
メッシュデータの人口総数を地図上に可視するには,シェープファイルと結合する必要がある.そこで,メッシュデータの境界データをダウンロードし,結合する.
統計地理情報システム > 境界データダウンロード > 4次メッシュ(500mメッシュ) > 世界測地系緯度経度・Shapefile > 都道府県で絞込みはコチラ > 13 東京都 > M5339
read.sf()を用いてデータを読込.meshと名付けるmeshとTokyo_popをKEY_CODEで結合.mesh<-
read_sf("HDDSWH5339/MESH05339.shp")
#データ結合
Tokyo_mesh<-
left_join(mesh, Tokyo_pop, by=c("KEY_CODE"))
#人口総数(T001101001)を数値(integer:整数)に変更
Tokyo_mesh %>%
mutate(T001101001=as.integer(T001101001)) ->
Tokyo_mesh地図上に可視化
Tokyo_meshは東京都全体の地図.ここでは,多摩市周辺だけを示したい.そこで,多摩市周辺の行政区域(Tama_shi_map)とTokyo_meshの重複する部分を抽出する.
st_intersection()を用いて,2つのシェープファイル(Tokyo_meshとTama_shi_map)の重複部分を抽出.
#準備
proj<-
"+proj=longlat +datum=WGS84"
Tokyo_mesh_proj<-
Tokyo_mesh %>%
st_transform(proj)
Tama_shi_proj<-
Tama_shi_map %>%
st_transform(proj)
st_agr(Tokyo_mesh_proj)="constant"
st_agr(Tama_shi_proj)="constant"
#多摩周辺の地図に人口が加わったメッシュデータの完成
Tama_mesh<-
st_intersection(x=Tokyo_mesh_proj, y=Tama_shi_proj)多摩ニュータウンと多摩周辺の人口分布の可視化
ggplot()
NA(欠損値)は情報(ここでは人口総数)が含まれないことを意味する.#完成図
ggplot()+
geom_sf(data=Tama_mesh, aes(fill=T001101001))+
scale_fill_viridis_c(option="G", direction=-1)+
geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
labs(fill="人口総数",
caption="出典:国土交通省国土数値情報")+
ggtitle("多摩NT位置と多摩市周辺人口(2020年)")+
theme_bw()mapview()
mapview()を用いる利点:マウスカーソルを自ら動かすことで知りたいメッシュ内の人口総数を調べられるようになる.3
layer.nameを用いて凡例を変更.#viridisカラーの使用(好み)
library(viridis)
#色の指示(好み)
my_color<-
viridis(7, option="mako", direction=-1)
#mapview表記
#col.regions=my_colorは好み(削除しても問題ない)
mapview(Tama_mesh, zcol="T001101001",
alpha.regions=0.6, col.regions=my_color,
layer.name="人口総数")+
mapview(Tama_newtown, cex=2.5,
color="yellow", col.regions="yellow")st_buffer():円(ポリゴン)を作成.ここでは,Tama_newtownのポイント(geometry列に納められている)を中心として500mの円を作成.bufferと名付ける.#単位を付与(ここでは周辺500m)
library(units)
Tama_newtown %>%
st_buffer(geometry, dist=set_units(500, m)) ->
buffer
#可視化1(500mの円)
ggplot()+
geom_sf(data=Tama_mesh, aes(fill=T001101001))+
scale_fill_viridis_c(option="G", direction=-1)+
geom_sf(data=buffer, color="red", fill="NA")+
geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
labs(fill="人口総数",
caption="出典:国土交通省国土数値情報")+
ggtitle("多摩ニュータウン周辺500mの範囲")+
theme_bw()st_intersects():Tama_meshから円(buffer)ポリゴンと重なるメッシュを検索.lengths()>0:
重なるメッシュを抽出しTRUEに.Buffer_Pop<-
lengths(st_intersects(Tama_mesh, buffer))>0 Tama_meshからBuffer_Pop==TRUEを満たすメッシュを抽出(filter()).
.ですべての変数を残す.Tama_mesh2と名付ける.Tama_mesh %>%
filter(., Buffer_Pop==TRUE) ->
Tama_mesh2#可視化2(円を重複するメッシュ(赤)に変更)
ggplot()+
geom_sf(data=Tama_mesh, aes(fill=T001101001))+
scale_fill_viridis_c(option="G", direction=-1)+
geom_sf(data=Tama_mesh2, color="red", fill="NA")+
geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
labs(fill="人口総数",
caption="出典:国土交通省国土数値情報")+
ggtitle("多摩ニュータウン周辺メッシュ")+
theme_bw()周辺500mのメッシュに含まれる人口総数の計算
distinct()で削除.
NA)が含まれる場合も統計値を計算できない.欠損値を無視して計算する(na.rm=TRUE).#重複するメッシュ(KEY_CODE)を取り除く
#(.keep_all=TRUE)その他の変数残す
Tama_mesh2 %>%
distinct(KEY_CODE,.keep_all=TRUE)->
Tama_mesh3
#欠損値無視オプション(na.rm)
Tama_mesh3 %>%
summarize(sum_pop=sum(T001101001, na.rm=TRUE))| sum_pop | geometry |
|---|---|
| 147078 | MULTIPOLYGON (((139.375 35…. |
#上記より下記のよう示す方がよいかも
sum(Tama_mesh3$T001101001, na.rm=TRUE)## [1] 147078
mapview()で確認
mapview()を用いて,周辺500mと重なるメッシュのみを可視化.
lwd:枠の太さを指示.legend)とホームボタン(homebutton)を削除(FALSE).#mapview
mapview(buffer, color="red",
alpha.regions=1, lwd=1.2,
legend=FALSE, homebutton=FALSE)+
mapview(Tama_mesh2, zcol="T001101001",
alpha.regions=0.6, col.regions=my_color,
layer.name="人口総数", homebutton=FALSE)+
mapview(Tama_newtown, cex=2.5,
color="yellow", col.regions="yellow",
legend=FALSE, homebutton=FALSE)謝辞
ポリゴンと重なるメッシュの検索,抽出方法について山本雅資氏にRコードを教示いただいた.
Rによる地理空間データの可視化
統計地理情報システムの国勢調査からは1995年から2020年(5年ごと)までの人口総数(男女),世帯総数などがわかる(アクセス日:2022-10-29).2015年からは年齢別人口,外国人人口など詳しく記録されている.↩︎
準備に示されているコードはコミュニティバス路線上の人口密度の可視化に説明あり.↩︎
mapview()を用いた地理空間データの可視化についてはインタラクティブな地図によるデータの可視化参照.↩︎